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The classical theory of electrokinetic phenomena is based on the mean-field approximation, that 
the electric field acting on an individual ion is self-consistently determined by the local mean charge 
density. This paper considers situations, such as concentrated electrolytes, multivalent electrolytes, 
or solvent-free ionic liquids, where the mean-field approximation breaks down. A fourth-order 
modified Poisson equation is developed that captures the essential features in a simple continuum 
framework. The model is derived as a gradient approximation for non-local electrostatics of interact- 
ing effective charges, where the permittivity becomes a differential operator, scaled by a correlation 
length. The theory is able to capture subtle aspects of molecular simulations and allows for simple 
calculations of electrokinetic flows in correlated ionic fluids, for the first time. Charge-density os- 
cillations tend to reduce electro-osmotic flow and streaming current, and over-screening of surface 
charge can lead to flow reversal. These effects also help to explain the suppression of induced-charge 
electrokinetic phenomena at high salt concentations. 



I. INTRODUCTION 

The classical theory of the electric double layer and 
electrokinetic flow near a charged surface is over a cen- 
tury old and remains in wide use today [![. The clas- 
sical theory has been extremely powerful in a number 
of diverse fields such as colloidal science, biophysics, mi- 
cro/nanofluidics and electrochemistry. While the useful- 
ness of the classical electrokinetic theory is not in ques- 
tion, there is a long history of recognizing the limitations 
and offering new formulations 0, y] • 

The equations are built on a set of assumptions which 
are clearly violated in various instances. The classical 
theory was developed for a surface in chemical equilib- 
rium with a dilute solution of point ions with a double- 
layer voltage on the order of the thermal voltage, kT/e — 
25 mV 0t6|. Stern recognized in 1924 that the assump- 
tion of point ions leads to predicted concentrations that 
are impossibly high at modest voltages. Stern introduced 
the idea of a molecular layer of finite size to reduce (but 
not eliminate) this un-physical divergence by imposing a 
distance of closest approach of ions to the surface [7[ . In 
many practical situations when the surface is unknown or 
uncontrolled, the macro-scale observable quantities such 
as capacitance or fluid slip velocity are fit with effec- 
tive Stern layer properties to bring the classic model into 
agreement with experiment. 

There has been recent interest in including finite ion 
size effects into the continuum electrokinetic model to go 
beyond the simple Stern layer approach Q. It is appar- 
ently not well-known that Stern proposed such an ap- 
proach as the final (un-derived) equation in his 1924 pa- 
per Q. One driver for interest in steric effects are appli- 
cations where electrokinetic phenomena are exploited in 
devices with electrodes placed in direct contact with the 
fluid f§-fllj|. These "induced-charge electrokinetic phe- 



nomena" [12| have shifted attention to a regime where 
double-layer voltage reaches several Volts « lOOfcT/e, a 
regime where the point ion theory is certainly invalid. 
To account for finite sized ions, a variety of "modified 
Poisson-Boltzmann equations" (MPB) have been pro- 
posed 0, [lH . The simplest possible MPB model is the 
one proposed (and subsequently forgotten) by Bikerman 
in 1942 jl4j . which is a continuum approximation of the 
entropy of ions on a lattice . Such modifications to the 
continuum theory can predict an otherwise unexplained 
high frequency flow reversal in AC electroosmotic pumps 
[16j . and capacitance of surfaces with no adsorption [2j. 

Extensions of the classical electrokinetic theory are 
also required for room-temperature ionic liquids (RTILs) . 
RTILs typically have large organic cations and similar or- 
ganic or smaller inorganic anions and hold promise as 
solvent-free electrolytes for super-capacitors, batteries, 
solar cells, and electro-actuators jl7H24j. For these appli- 
cations, data for the RTIL/metal interface has typically 
been interpreted through models based on the classical 
theory despite the fact that this dense mixture of large 
ions bears little resemblance to a dilute solution of point- 
like ions. Recently, Kornyshev [25[ stressed the impor- 
tance of finite-sized ions and developed a theory equiva- 
lent to Bikerman's, where the bulk volume fraction can 
be tuned to describe electrostriction of the double layer. 

In spite of some success in applying a theory which ac- 
counts for steric hinderance in electrolytes at high voltage 
and RTILs, these models are unable to describe short- 
range Coulomb correlations [26| . In many important 
situations, classical theory breaks down due to strong 
correlations between nearby ions. In concentrated so- 
lutions, systems with multivalent ions (relevant for bi- 
ology), RTILs, or molten salts, electrostatic correlations 
which go beyond the mean electrostatic potential become 
dominant. Correlations generally lead to over-screening 
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of a charged surface, where the first layer provides more 
counter-charge than required; the next layer then sees a 
smaller net charge of the opposite sign, which it over- 
screens with excess co-ions; and so-on. 

Such overscreening is usually studied with molecular 
dynamics simulations, Monte-Carlo simulations (MC), 
Density Functional Theory (DFT), or integral equation 
methods based on the statistical mechanics of charged 
hard spheres. While these simulations are based on more 
realistic assumptions than classic theory, the complexity 
prohibits analytical progress and the computational cost 
and complexity can be high. In many applications we 
are interested in charging dynamics, fluid flow, or other 
macroscale behavior where a simple model is needed. To 
date, essentially all modeling of electrokinetic flow has 
been based on the mean field approximation, where the 
electric field acting on the ions is self consistently deter- 
mined by the mean charge density. 

In this paper we maintain a continuum formulation 
and develop a modified Poisson equation which accounts 
for electrostatic correlation effects in diffuse electric dou- 
ble layers. This model is applicable to concentrated or 
multivalent electrolytes, room temperature ionic liquids, 
and molten salts. Recently, we (along with A. A. Ko- 
rnyshev) derived and applied this continuum model for 
RTILs 27]. In that work, we found good agreement in 
terms of the double layer structure and the capacitance 
when compared to molecular dynamics simulations. In 
the present work, we present the derivation in detail and 
apply the same continuum model to electrolytes, where 
correlations become important at high salt concentration 
and with multi-valent ions. We also couple the modified 
electrostatic theory to the Navier-Stokes equations, as 
we (along with M. S. Kilic and A. Ajdari) recently pro- 
posed 0. From this theoretical framework, we compute 
electrokinetic flows beyond the mean-field approximation 
for the first time. The model predictions are also com- 
pared to molecular simulations and some experimental 
data. 

Before we begin, we emphasize that any attempt to de- 
velop and modify continuum models for molecular scale 
phenomena is fundamentally limited. Nevertheless, our 
goal is to develop and test models that are simple enough 
to facilitate a better understanding of electrokinetics in 
macroscale experiments and devices. In particular, we 
describe flows in correlated electrolytes and ionic liquids 
with only one new parameter, an electrostatic correlation 
length. 



II. CONTINUUM ELECTROKINETIC 
EQUATIONS 

A. Classical mean-field theory 

The classic theory of electrokinetics assumes a dilute 
solution of point ions. The electrochemical potential, /Zj, 



of the i ionic species in an ideal dilute solution is, 

^idcal = kThgCi + z . e(j) (1) 

where k is Boltzmann's constant, T is the temperature, 
Ci is the concentration, zi is the charge number, e is the 
elementary charge and <fi is the electric potential. We 
relate the flux of each species, F i; to the gradient in the 
chemical potential and conservation of mass yields, 

_2=_ V .F i = -V- Uu-^ Ci VMi,J. (2) 

where Di is the diffusivity and u is the mass averaged ve- 
locity. It is important to remember that directly relating 
the flux of each species to its own gradient in chemical 
potential is an assumption that is strictly only valid in 
dilute solutions. This relationship assumes that the dif- 
fusivity tensor is diagonal. The system is traditionally 
closed by making the mean field approximation in which 
the electric potential satisfies the Poisson equation, 

- V • eV(j) = p = 2J ZitCi, (3) 

i 

where p is the charge density and e is the permittivity. 
Equations are typically referred to as the Poisson 
Nernst Planck (PNP) equations. The PNP equations are 
coupled to the Navier Stokes (NS) equations for fluid 
flow, where an electrostatic force density, pV0, is added, 

p m (j£ + u ■ Vuj = -Vp + r/V 2 u - pV0, (4) 

Vu = 0, (5) 

where r\ is the viscosity, p m is the mass density, and P is 
the pressure. In the classical theory the fluid properties 
such as the viscosity and permittivity are usually taken 
as constants. 

Solutions to equations [5] - E] require boundary condi- 
tions. Boundary conditions can vary depending on the 
physical situation. Typically, the no-slip condition for 
fluid velocity is assumed, but modifications can allow for 
slip at a solid surface. A common boundary condition 
for the ion conservation equation is that there is no flux 
of ions at a solid surface. However, in cases with elec- 
trochemical reactions or ion adsorption, other boundary 
conditions are required. 

The boundary condition for the potential depends 
upon the physics of the interface. Our interest is on metal 
electrode surfaces where one can simply fix the applied 
potential 4> = 0o or allow for a thin dielectric layer (or 
compact Stern layer) on the electrode surface through 
the mixed boundary condition (28| . 

A0s = 4> ~ 0o = ^sn ■ V0 - (6) 

where As = ehs/es is an effective thickness of the layer, 
equal to the true thickness h$ multiplied by the ratio 
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of permittivities of the solution e and the layer es, and 
Cs = £ s/h,s is its capacitance. When applying Q to a 
metal electrode, one can set q$ — to model the Stern 
layer as a thin dielectric coating of solvent molecules [2!| , 
while specific adsorption of ions would lead to qs 7^ 0. 

While the PNP+NS formulation is widely studied and 
widely used, the mathematical solution can be compli- 
cated. In many cases we can make mathematical sim- 
plifications that allow for analytical progress or simple 
models to be derived from the PNP+NS starting point. 
In this work, we are developing a physical modification 
to the equations. 



B. Modifications for chemical effects 

In a recent review article we (along with M.S. Kilic and 
A. Ajdari) discuss in detail a number of ways in which 
the classical mean-field theory of electrokinetics breaks 
down and propose some simple modifications for large 
voltages and concentrated solutions We stress that 
attempts to go beyond the classic equations have a long 
history and refer the interested reader to Ref. |2] for a 
more complete account of the literature. 

To account for various thermodynamic non-idealities 
in concentrated solutions, we can extend the chemical 
potential by adding an excess term to that of the ideal 
solution, 

Hi = kTlogCi +Zie(j) + [i ex . 

In the case of volume constraints for finite-sized ions, 
following Bikerman this excess chemical potential 
could be written as 



-fcTln(l - $) 



(7) 



where $ is the local volume fraction of ions. The same 
model of the excess chemical potential can also be derived 
from the configurational entropy of ions in a lattice gas 
in the continuum limit, as first noted by Grimley and 
Mott [15].. We attribute this model to J.J. Bikerman 
though it has been independently rediscovered at least 
seven times since then and was possibly first discussed by 
Stern in 1924. Other approaches can be used to modify 
the chemical potential for volume constraints, such as 
Carnahan Starling equation of state for the entropy of 
hard-spheres in the local density approximation [30l432l |. 
Regardless of the model for steric volume constraints, 
these modifications all allow the formation of a condensed 
layer of ions very close to the surface at high voltage. 
This layer forms at high voltage as the classic theory 
allows for an impossibly high density of ions. 

Another modification we have discussed in detail is 
charge induced thickening, where one supposes that the 
viscosity of the fluid depends upon the local charge den- 
sity and typically increases in the inner part of the double 
layer (effectively moving the "shear plane" of no slip away 
from the surface). Charge induced thickening provides 



a possible explanation for the decay in induced-charge 
electro-osmotic flow 12j that is observed in many ex- 
periments at high salt concentration and/or high voltage 
js| H3[. Below, we will argue that electrostatic correla- 
tions may also play a significant role in explaining the 
data. 

The permittivity e of a polar solvent like water is usu- 
ally taken as a constant in ([3]) , but numerous models exist 
for field-dependent permittivity e(|V0|), as discussed in 
Q. The classical effect of dielectric saturation reduces 
the permittivit y a t large fields due to the alignment of 
solvent dipoles [29|, |3J-[36( , although an increase in dipole 
density near a surface may have the opposite effect [37j . 
A recent model which included excess ion polarizability 
demonstrated excellent agreement with experimental ca- 
pacitance data on surfaces with no adsorption [381 ] . 

While these and many other modifications have been 
explored, in this work we only consider the additional 
chemical effect of finite ion size, so we can focus on novel 
effects of electrostatic correlations. 



C. Simple modification for Coulomb correlations 

The most fundamental modification of the classical 
theory, which has resisted a simple treatment, would be 
to relax the mean-field approximation. While the study 
of electrostatic correlations in electrolytes has a long his- 
tory, we are not aware of any attempts to go beyond 
the mean-field approximation ([3]) in dynamical problems 
of ion transport or electrokinetics. Dynamical problems 
with bulk flow would seem to require a simple contin- 
uum treatment of correlation effects, ideally leading to a 
general modification of Eq. |3J 

In recent work on RTIL, we (along with A. A. Korny- 
shev) derived a Landau-Ginzburg type continuum model 
which accounts for electrostatic correlations in a very 
simple and intuitive way [27]. A general derivation based 
on nonlocal electrostatics will be developed in the next 
section, but first we present the final result, which is a 
modified fourth-order Poisson equation, 
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and a modified electrostatic boundary condition, 
n D = n ■ e(^V 2 - 1)V<£ = q s 



(8) 



(9) 



where D is the displacement field. Due to Coulomb cor- 
relations, the effective permittivity i, defined by D = 
— eV0, is a linear differential operator, 



s = e (l - £ 2 V 2 



(10) 



This unusual dielectric response, signifying strong corre- 
lations, is consistent with some well known properties for 
molten salts, although we extend it here to more general 
situations. In particular, for small, sinusoidal perturba- 
tions of the electric field of wavenumber k, the corre- 
sponding small- k expansion of the Fourier transform of 



4 



the permittivity, 

e(k) = e [1 + (4fc) 2 ] (11) 

grows with k in the case ao > where correlations 
promote charge density oscillations and discrete cation- 
anion-cation-anion... ordering. This matches known re- 
sults for molten salts, although at smaller wavelengths 
(larger k) the permittivity transform e(k) has diver- 
gences due to electronic relaxation and other phenom- 
ena [39|, Here, we do not use the notion of 
wavelength-dependent permittivity, which only applies to 
small periodic bulk perturbations. Instead, we introduce 
the concept of a permittivity operator in Poisson's equa- 
tion, which can be applied to general nonlinear response 
in asymmetric geometries and near surfaces. The new 
parameter £ c is an effective length scale over which cor- 
relation effects are important, discussed below. Its value 
is not precisely known, though we can place approximate 
bounds on its value. 

Similar higher-order Poisson equations have been de- 
rived as approximations for the equilibrium statisti- 
cal mechanics of point-like counterions (one component 
plasma) near a charged wall [4l| - |43l |; Santangelo (4l| 
showed that © is exact for both weak and strong cou- 
pling and a good approximation at intermediate coupling 
with £ c set to the Bjerrum length; Hatlo and Lue [43[ de- 
veloped an approximation for £ c . The extension to elec- 
trolytes and non-ideal solutions was first proposed in our 
review paper Q as part of a general modeling framework 
for electrokinetics, but without a derivation or any ex- 
ample calculations. In our recent work on RTILs (27J, we 
presented a general variational derivation of the model 
and first applied this modified Poisson equation to pre- 
dict double layer structure and capacitance (RTIL), using 
the ion size as the correlation length scale. 

Since Poisson's equation © is now fourth-order, we 
need an additional boundary condition. For consistency 
with our derivation below, we neglect correlations very 
close to the surface (at the molecular scale) and apply the 
standard boundary condition, —en ■ V0 = q s . Equation 
([9]) then implies 

n-V(V 2 ?y)=0 (12) 

which requires that the mean-field charge density is "flat" 
at the surface. Although this boundary condition is 
consistent with the derivation of our model, we stress 
that it is neither unique nor rigorously derived. Alter- 
native boundary conditions should be considered in the 
future, including the possibility of nonlocal models (e.g. 
which are required to describe density oscillations result- 
ing from packing constraints 26}). Here, we use Eq. ([12")) 
partly for its elegant simplicity and partly since it seems 
to consistently produce remarkably good results with our 
model in comparison to molecular simulations, not only 
for RTIL but also for concentrated electrolytes, as 
described below. 



III. DERIVATION OF THE MODIFIED 
POISSON EQUATION 

The following derivation is adapted from the support- 
ing information of our recent publication with A. A. Ko- 
rnyshev [27[, providing some more details and explana- 
tions of the steps. 



A. Electrostatic energy functional 

We begin by postulating general free energy functional 
broken into chemical and nonlocal electrostatic contribu- 
tions. Let G = G e i + Gchem, where G e i is the electro- 
static energy and G c hem = f v drg is the chemical (non- 
electrostatic) part of the total free energy, G. Suppose 
that Gchem is known, and let us focus on electrostatic 
correlation effects in G e i ■ 

The electrostatic potential, cf>, is defined as the elec- 
trostatic energy per ion (free charge). The electrostatic 
energy cost for adding a charge bp in the bulk liquid vol- 
ume V or 5q s on the metal surface S is, 

SGei = [ drcj)Sp+ [ drcj)Sq s . (13) 
Jv Js 

The charge is related to the displacement field D via 
Maxwell's equation, 

V ■ D = p => Sp = V -SB. (14) 

The corresponding boundary condition for an ideal metal 
surface (where D = 0) is, 

[n - D] = n ■ D = -q s => 5q s = -h ■ D. (15) 

Substituting these expressions into Q13[) and using Gauss' 
theorem, along with the definition of the electric field, 
E = — V0, we recover the standard electrostatic free en- 
ergy equation |44| . 

SGei = [ drE-SB. (16) 
Jv 

In the linear response regime (for small external elec- 
tric fields), we have 

D = eE, (17) 

where e is a linear operator, whose Fourier transform 
e(k) encodes how the permittivity depends on the wave- 
length 2n/k of the k- Fourier component of the field, due 
to discrete ion-ion correlations, as well as any non-local 
dielectric response of the solvent. A crucial feature of our 
approach, however, is that we do not restrict ourselves to 
small amplitude perturbations in Fourier space. Instead, 
we consider a general linear permittivity operator in real 
space and focus on correlation effects. 
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By linearity, we can integrate (|16l) over <5D through a 
charging process that creates all the charges in the bulk 
and surface from zero to obtain 



G e 



drED. 



(18) 



For a given distribution of charges p and q s , with asso- 
ciated displacement field D, the physical electric field 
E is the one that minimizes G e i, subject to the con- 
straint of satisfying Maxwell's equations (|14[) - (|T5|) . Since 
E = — V</> to enforce V x E = 0, we can minimize G e i 
with respect to variations in cf>, using Lagrange multipli- 
ers Ai and A2 to enforce the constraints, 



^E • D + Ai (p- V • D) 



G e/ [0] = f dr 
Jv 

+ <j> dr s A 2 {q s + h ■ D) 



(19) 



To calculate the extremum, we use the Frechet functional 
derivative: 

SGa = Um G el [cl) + e<f)o8 e }-G e i[(f)} 
5(f> e->-o e</>o 

where S<fi e = o <5 e (r, r') is a localized perturbation of the 
potential (with compact support), which tends either to 
a 3D delta function in the liquid (r € V) or to a 2D delta 
function on the surface (re S) as e — > 0, and 0o is an 
arbitrary potential scale for dimensional consistency. By 
setting 8G e i/8<j) = for both surface and bulk variations, 
we find Ai = A2 = 4>- Finally, using vector identities, we 
arrive at a general functional for the electrostatic energy, 



G d [(t>] = J^dr (pcj) + • + j> dr s q t 



4, (21) 



whose variational derivative with respect to 4> will be set 
to zero, once we know the relationship between D and 
E = -V<t>. 



B. 



Nonlocal electrostatics for correlations 



To model the field energy, we assume linear dielectric 
response of the individual molecules (ions and solvent) 
with constant mean permittivity e, plus a simple non- 
local contribution for Coulomb correlations. Here, the 
permittivity e describes the electronic polarizability of 
the ions (for RTIL) as well as (in the case of electrolytes) 
the dielectric relaxation of the solvent. There is an ex- 
tensive literature on nonlocal electrostatic models of the 
form, D(r) = J dr'e(r, r')E(r'), mainly focused on de- 
scribing the nanoscale dielectric response of water (45l — 
l48j . In this work, we take a very different approach 
because our aim is to model the transient formation of 
correlated ion pairs of opposite sign (zwitterions) , which 
act as dipoles and contribute to the nanoscale dielectric 
response of strongly correlated ionic liquids. 



The theory begins by postulating a non-traditional 
form of the energy density stored in the electric field in 
the dielectric medium, 
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-~V<£-D = | (E(r) 2 + J dv'K{v, v')p{r)p{v' 

(22) 

where we define 



eV ■ E 



-eV 2 



(23) 



as the "mean-field charge density" , which would produce 
same the electric field in the dielectric medium without 
accounting for ion-ion correlations |27j ]. In this theory, 
nonlocal electrostatic effects are assumed to derive from 
pairwise interactions between effective charges, defined 
in terms of the local divergence of the electric field via 
the standard second-order Poisson equation with con- 
stant permittivity e. The non-local kernel K(r,r') is in- 
tended to describe correlations between discrete pairs of 
fluctuating ions resulting from Coulomb interactions in 
the liquid. 

To take into account screening phenomena, we assume 
that the correlation kernel K (r, r') decays exponentially 
over a characteristic length scale £ c . Below this distance, 
ions experience bare Coulomb interactions, and beyond 
it, thermal agitation and many-body interactions act to 
suppress direct electrostatic correlations. The correlation 
length is clearly bounded below by the ion size a, which 
becomes the most relevant length scale in a highly con- 
centrated electrolyte (including the solvent shell in the 
ion size) or a solvent-free ionic liquid. In the simplest 
version of our theory for dense ionic mixtures, it is possi- 
ble to avoid adding any new parameter by simply setting 
£ c = a, as in our work on RTIL [27| ■ In dilute elec- 
trolytes, however, the correlation length should increase 
with the mean ion spacing, and we expect it to be cut 
off at the scale of the Bjerrum length £b, which is the 
separation distance between ions below which the bare 
Coulomb energy exceeds the thermal energy ksT. 

In order to obtain a simple continuum model, we fur- 
ther assume that charge variations mainly occur over 
length scales larger than £ c (corresponding to small per- 
turbation wavenumbers, £ c \k\ <C 1). In this limit, we 
perform a gradient expansion for the non-local term 
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|V0| 2 



n=0 



£T 



-V n p 



(24) 



where a n are dimensionless coefficients, which implies 



G, 



\V4>\ 2 + Y j a n . 2 {£ n c - x V n (l 



n=2 



dr s q s 



(25) 



For simplicity, we typically truncate the expansion after 
the first term, even though this may become inaccurate 
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in situations of interest with charge density variations at 
the correlation length scale. 

From the gradient expansion of the nonlocal electro- 
static energy functional, we set SG e i/5(j) = for bulk and 
surface perturbations in (|25[) . In this way, we recover 
Maxwell's equations (fTi |) -([T5 |l . with 



(26) 



D = eE, 



where the permittivity operator has the following gradi- 
ent expansion, 



, a n -i 



(27) 



Eq. (|10p results from the first term in the gradient ex- 
pansion with the choice ao = 1 (after suitably rescaling 
i c ), where the overall negative sign of this term is chosen 
to promote over-screening. The corresponding small-fc 
expansion of the Fourier transform of the permittivity, 



1 + ^a n _i(- 



l (4fc) 2 



(28) 



grows with k at small wavenumbers in the case where cor- 
relations promote overscreening, «o > 0, as noted above. 
This is a hallmark of Coulomb correlations, promoting 
alternating charge ordering. 



IV. CORRELATED ELECTROKINETICS AT A 
PLANAR SURFACE 

A. Basic equations 

To demonstrate how correlation effects may influence 
double layer structure and electrokinetic flows, we start 
by exploring the behavior at a planar surface. We assume 
a ID double layer at equilibrium with constant e and a 
z + : z~ binary electrolyte. 

The model we solve is thus, 



To calculate hydrodynamic slip, we start with the 
Navier Stokes equation and assume that in the electric 
double layer we have a balance between the electric body 
force and viscous forces, 



= 7] 



cPu 
dx 2 



PeE t 



where E t is the electric field tangential to the surface. In 
our model this becomes, 







d 2 u 



dx 2 



Et 



As with the standard Hclmholtz-Smoluchowski equation, 
we can integrate this equation across the double layer 
twice to obtain (with the convention that far from the 
wall, 0=0). 



u(oo) = 



T] 



1 



il d 2 



0(0) dx 2 



x=0 



In the above expression, we are assuming that the 
medium permittivity and viscosity are constant within 
the double layer, though this approximation can be re- 
laxed. An important general prediction is that the 
classical Helmholtz Smolukowski slip velocity, Uhs = 
—eE t <j)(0)/rj, is modified by the inclusion of correlation 
effects. This can be understood as a consequence of 
nonuniform permittivity. 

The total charge in the double layer is given as the 
integral of the charge over the double layer, 



Q : 



p e dx = I e 



c dx± dx 2 



Evaluating this integral and using the boundary condi- 
tions at a solid electrode stated above we obtain, 



p e dx = e — — 
dx 



x=0 



with the total capacitance defined as C = q/<p(Q). 



: dx* 



dx 2 



= p = z 



+ 

ec — z 



The boundary conditions at the electrode surface of fixed 
potential are, 

This electric potential equation is solved along with the 
equations that the chemical potentials must be constant 
at equilibrium. In this work we consider the Biker- 
man model for volume constraints only with equal sized 
cations and anions such that the chemical potential of 
the ions is, 

[X± = fcTlogQ - fcTlog(l - a 3 (c + + c_)) ± z±e<f) 



B. Dimensionless formulation 

We assume a binary z + : z~ electrolyte such that the 
far field concentrations of the cations and anions follows 
£ + c+ = z^c^. For simplicity we assume that the cations 
and anions are of the same diameter. We make the for- 
mulation dimensionless using the scales c + = c + /c+, 
c~ = c^/c+ , and 4> = 4>{e/kT). The dimensionless con- 
centrations can be written as explicit functions of the 
electric potential, 



Z (j> 



(29) 



(30) 
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where the function, /3, is given by 



i 



z-+z+ 



z - e -z + 4> _|_ z + e z <f> 



(31) 



where v is the volume fraction in the bulk and has a value 
v = (1 + §^)c+ a 3 . For the case of a 1:1 electrolyte note 

that P{4>) = l/(l + z/(cosh(>)-l)) = l/(l + i/sinh 2 (^/2)) 
as has been used in previous works [Hi, .25]. We relate the 
lattice size parameter, a, to the spherical ion diameter, 
d, as a 3 = §d 3 /0.63 = 0.83d 3 where the factor of 0.63 is 
the maximum volume fraction for random close packing 
of spheres. 

The Poisson equation is scaled by the Debye length; 
i.e. x = x/\d where 



Xd = 



ekT 



e 2 ct c z+(z+ + z ) 
Under this scaling our governing equation becomes, 



d 2 cf> ? 2 d A 4> 
dx 2 ' "di 1 



_ g-Z + (/> _ gZ (f> 



(z+ + z-) 



(32) 



where 6 C = £ c /Xd- This equation is subject to the bound- 
ary conditions that the potential at the electrode is fixed, 
the third derivative of the potential is zero, and the po- 
tential goes to zero smoothly at infinity. 

There are three dimensionless parameters which 
emerge from our formulation, the bulk volume fraction 
v, the correlation length scale 8 C , and the applied volt- 
age (or known surface charge) . The solution also depends 
on the valences of the ions z + and z~ . 

In dimensionless terms, the slip velocity relative to the 
Helmholtz-Smulokowski velocity is, 



u{oo) 
Uh7 



8 2 d?4> 



c/>(0) dx 2 



(33) 



where Uhs — sEt<p{0)/rj. The capacitance relative to the 
Debye-Huckel capacitance, Cdh = e/Ad is simply, 



C 



C 



DH 



0(0) dx 



(34) 



5=0 



For the remainder of the paper we will drop the tilde 
notation and only refer to dimensionless quantities in our 
equations. 



C. Low voltage analytical solutions 

When the voltage is small relative to the thermal volt- 
age, kT/e, the problem is drastically simplified and the 
right hand side of our equation becomes, 



d 2 ^ 
dx 2 



c dx± 



(35) 



This equation can be solved analytically, though the form 
depends upon whether the value of S c is less than, greater 
than, or equal to 1/2. 



1. Solution 8c < 1/2, "weak" correlation effects 

When S c < 1/2 the analytical solution at low voltage 
has the form, 



4>{x) 



i-ki/k 



h 3 

2 



(36) 



where 



fcl = V — W 2 — ' fc2 = V — — ■ 

The capacitance of the double layer is, 

c MHf) 
c DH i-kim ' 

and the slip velocity is, 



(37) 



Uhs 



l-k\/k. 



(38) 



In the limit of 8 C going to zero k\ = 1 and ki = oo, thus 
we recover the Debye-Huckel solution <f>{x) = <fi(0)e~ x . 
This new solution has a functional form very similar to 
the classic double layer. The structure is given as the 
sum of two exponentials with decay lengths on the order 
of unity, though slightly modified. 



2. Solution for S c > 1/2, "strong" correlation effects 

When S c > 1/2 the analytical solution at low voltage 
has the form, 

4>{x) = ct){Q)e- klX (cos(k 2 x) - A sm(k 2 x)) (39) 

where 



VWTT ^28 c - 1 V25 C + 1(<5 C - 1) 

fc l = 7Tc > fc 2 = TTc j A = 



2S C 28c 1 V25 c - l(<y c + l) 

The capacitance of the double layer is, 



C y/28 c + 1 



Cdh S c + 1 



(40) 



which decays with increasing correlations. The slip ve- 
locity is, 



Uhs V S c + 1 



(41) 
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a) b) c) 




FIG. 1: Low voltage solutions to the continuum model, a) Double layer structure at different values of S c . Solutions are shown 
for S c = (dashed) and S c — 1, 2, and 5. b) Capacitance and c) slip velocity as a function of the correlation length scale, S c . 



The slip velocity changes sign if S c is sufficiently large. 
In particular, there is an electro- osmotic flow reversal or 
electrokinetic charge inversion of the surface when the di- 
mensionless correlation length exceeds the golden mean: 
S c > (1 + a/5)/2 = 1.618. 

The form of the double layer becomes modified as S c 
increases. We find that the functional form consists of 
decaying sinusoids with a length scale provided explicitly 
by fci and ki- At relatively large values of S c the length 
scale of the decay and the oscillations is approximately 



D. Numerical results 

At low voltage, the solution has only one free parame- 
ter, the correlation length scale, S c . The structure of the 
double layer as S c is varied is shown in Figure [1] We see 
that as the strength of the correlations is increased the 
double layer shows charge density oscillations. From the 
analytical solution we see that the oscillations emerge 
when 5 C is greater than 1/2. The length scale for the 
whole double layer also increases as the correlations are 
increased. From the analytical solution we can easily see 
at large 8 C that the size of the double layer grows with 
the square root of S c . For small values of 8 C the results 
become indistinguishable from the classic Debye-Huckel 
solution. 

In Figure [1] (b) we show the capacitance and (c) slip 
velocity as a function of Sc- We see a decrease in the 
slip velocity and the capacitance with increasing S c . As 
correlation effects become stronger the flow is quenched 
and then reverses direction. Note that from the analyti- 
cal solution that at S c = 1 that the flow velocity is half 
of Uhs and the flow reverses when S c > 0.618. These 
values of S c are easily reached at high concentration in 
aqueous electrolytes, as we will soon see. 

At higher applied voltage the structure of the solution 
changes dramatically as we show in Figure [5] Here we 
show sample solutions for a 1:1 and 2:1 electrolyte of 0.3 
nm ions as the voltage is changed. In Figure^ we show 



the structure of the double layer at different voltages at a 
cation concentration of 1 molar. Using the ion size as the 
correlation length scale and as the volume fraction then 
for the 1:1 system S c = 0.988, v = 0.0270 and for the 2:1 
system S c — 1.71, v = 0.0405. As the voltage increases, 
the charge density at the wall saturates to a value deter- 
mined by the steric constraints. This condensed layer of 
ions grows as the voltage is increase. Without the correla- 
tions effect the charge density would decay monotonically 
from the maximum value to zero far from the wall. How- 
ever, with the correlation effects included in the model, 
the charge density oscillates and changes sign. These os- 
cillations are more pronounced in 2:1 system when the 
divalent ions crowd the wall. 

Turning to the capacitance in Figure [2h we find that 
correlation effects reduce the capacitance. The dimen- 
sionless capacitance is always 1 at zero voltage when 
5 C = 0, however when 5 C > the capacitance at zero 
voltage reduces according to Figure [TJd. At higher volt- 
age, the shape of the capacitance curve is similar to when 
S c = 0, the values are simply lower. This reduction in ca- 
pacitance is consistent with previous work on steric con- 
straints with the Bikerman model which found generally 
that the theory needed ion sizes that were bigger than 
one would expect physically to fit the experimental data 

urn 

The most dramatic departure from the classical model 
comes when computing the slip velocity in Figure [5J;. 
We see that at high concentration the model can predict 
reverse flow even at small voltages in the 2:1 system. At 
low concentration, we find that the model predicts classic 
slip at low voltage but predicts reverse flow as the voltage 
is increased even moderately. As the voltage is increased 
further, the model predicts the forward component of 
the flow begins to increase as the condensed layer grows. 
At high voltage the slip velocity for all concentrations 
begin to come together as the condensed layer begins to 
dominate the double layer structure. 

These preliminary flow results must be interpreted 
with caution. The model currently does not account for 
changes in the viscosity of the solution near the wall in 
the condensed layer. It is also unclear (as it is in classi- 
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a) b) c) 




x/a V V 

FIG. 2: Example solutions for a 1:1 (a, b, and c) and 2:1 electrolyte (d, e, and f) with 0.3 nm ion sizes, (a) Double layer 
structure showing the charge density profiles at wall voltages of -1, -2, -5 and -10 in units of kbT/e for a 1 Molar concentration 
of cations, (b) Dimensionless capacitance as a function of voltage for concentration of 0.01 (black), 0.1 (red) and 1 (blue) 
molar from top to bottom. Corresponding solutions with no correlation effects (5 C = 0) are shown as the dashed lines, (c) 
Dimensionless slip velocity as a function of voltage for concentration of 0.01 (black), 0.1 (red) and 1 (blue) molar. Without 
correlations the slip velocity is always 1. Figures (c), (d), and (e) are the same, only for the 2:1 system. The asymmetry is 
easily seen in the capacitance and slip velocity. 



cal theory) where the slip plane should be placed. Recent 
work by Jiang and Qiao shows via molecular dynamics 
simulations that electroosmotic flow can be amplified by 
short wavelength hydrodynamics 49]. These effects (and 
others) are not included in our model and may be re- 
quired for more detailed comparisons with experimental 
data. 



V. VALIDATION 
A. Comparison with molecular simulations 

In order to determine whether this model has validity 
in the context of aqueous electrolytes, we can compare 
the model predictions to those made by more sophisti- 
cated simulations such as Monte Carlo or Density Func- 
tional Theory (DFT). Monte Carlo simulations are often 
considered the standard for equilibrium chemistry while 
DFT has proven to quantitatively compare well against 
Monte Carlo at a much lower computational cost (Hoj . 
Our aim is to determine whether an even simpler contin- 
uum model can capture the same features. 

In a prior paper we compared this correlations model 
to molecular dynamics simulations of ionic liquids [27| . 



In that work we assumed that the potential at x = in 
the continuum theory was the potential offset from the 
wall by the radius of the ion. In comparisons to data 
for electrolyte solutions that follow, we find that here we 
obtain good results by taking the voltage at x — to be 
the electrode, i.e. not accounting for the radius of the 
ion as it approaches the surface. 

In Figure [3] we compare the ion distributions, g(x) — 
c(x) /coo , predicted by the continuum model to those pre- 
dicted by Monte Carlo simulations of Boda et al . The 
conditions here are a 2:1 electrolyte with surface charge 
of -0.3 C/m 2 and an ion diameter of 0.3 nm. We find that 
the continuum model predicts much of the same struc- 
ture as the Monte Carlo simulations, though the length 
scale of the oscillations and the amount of over- screening 
predicted by the continuum model is larger than seen in 
the simulations. Better agreement can be obtained by 
reducing the correlation length scale by about 50 per- 
cent. However, the classic electrokinetic model can only 
predict ion profiles which decay monotonically, so it is 
interesting that this extension for correlations effects can 
provide the basic double layer structure with no fitting 
parameters. 

In Figure U a-b we compare the continuum model to 
the Monte Carlo simulations looking at the relationship 
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FIG. 3: Comparison of the continuum model (solid lines) to 
Monte Carlo simulations of ref [5lJ . The conditions here are 
a 2:1 electrolyte with surface charge of -0.3 C/m 2 and an 
ion diameter of 0.3 nm. The points are the Monte Carlo 
simulation and the solid lines are the continuum model. 

between the double layer charge and electrode voltage. In 
a) we show results for a monovalent ion and in b) we show 
a 2:1 electrolyte at two different concentrations. The 
continuum model predicts the basic trends of the more 
complicated MC simulations, though under-predicts the 
voltage for a given charge. The inclusion of correlation 
effects brings the continuum results in better agreement 
with the MC simulations than when we only account for 
finite size effects. 

In Figures [5] we compare the continuum model to re- 
sults of density functional theory (DFT) simulations of 
Gillespie et al 50] for a 2:1 electrolyte. In Figures [5] we 
show curves of constant voltage over a range of surface 
charge and concentration. The results with the contin- 
uum model are in reasonable agreement with the DFT re- 
sults, especially at large concentrations and high charge. 
Importantly, the shape of these curves computed with 
DFT are well predicted by this simple continuum model. 
When S c = and at high concentration, the continuum 
model qualitatively departs from the DFT results. What 
is interesting about the continuum model with correla- 
tions included is that there are no fit parameters. 



B. Comparison with experiments 

We can also compare the model to an experiment, 
rather than to other simulations, as a more definitive 
test. In Figure|H]we compare the model to the nanochan- 
nel electrokinetic data of van der Heyden et al. [52j as 
was done by Gillespie et al In the experiment a 

nanochannel with a characterized surface charge is driven 
by a pressure driven flow and the streaming current is 
measured. In this case the flow is driven by pressure and 
not electro-osmotically. To compute the streaming cur- 




0.05 0.1 0.15 0.2 0.25 0.3 
Q (a 2 /e) 
b) 




-6 



-0.2 -0.1 0.1 0.2 0.3 

Q (a 2 /e) 

FIG. 4: Comparison of the continuum model with correlations 
(solid lines) to Monte Carlo simulations of ref 51] (points) and 
the continuum model with only steric effects (dashed lines). 
The ion diameter is 0.3 nm. In a) we show the result for a 1:1 
electrolyte and in b) we show the result for a 2:1 electrolyte. 
The upper solid blue curve and dots is for 0.1 M and the lower 
red curve and '*' is for 1 M concentration. 

rent we simply multiply our charge density profiles by 
the pressure driven velocity profile; 

I = W f p{x)u{x)dx (42) 
Jo 

where W is the channel width of 50 /xm, H is the chan- 
nel height of 450 nm, and u is the parabolic velocity 
profile. Since the double layer is so thin relative to the 
channel height of 450 nm, we can safely assume that the 
pressure driven velocity profile is locally linear at the 
wall; du/dx = AAPH/(Lri) for Pouiselle flow. Thus to 
compute the current per unit pressure drop for pressure 
driven flow we calculate, 

I 4WH f°° 
AP = —J ^> dx - ™ 
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FIG. 5: Comparison of the continuum model accounting for 
correlation effects (blue solid lines) to the DFT simulations 
of Gillespie et al [5(j] (black dashed) to the continuum model 
with 6 C = (red dotted lines). The ion diameter is 0.3 nm in 
the models and DFT. 




C(M) 

FIG. 6: Comparison of the continuum model experimental 
nanochannel data of [5^ ]. The electrolyte is a 2:1 with an 
assumed ion size of 0.3 nm. 



The current per unit pressure as a function of concen- 
tration is plotted in Figure |6] comparing the continuum 
model to the experiment. The agreement is qualitatively 
correct and predicts a reversal in the current around the 
same concentration as seen in the experiments. The 
slower velocities at high concentration seen in the ex- 
periment is consistent with charge induced thickening, 
and increase in viscosity in a condensed layer of ions [2| . 
There is still uncertainty in application of this model for 
flow. It is unclear where the slip plane should sit and 
whether the solution viscosity near the wall should be 
taken as a constant. This uncertainty applies equally to 
the continuum model and the DFT simulations, as in 
those simulations the current is calculated in the same 
way it is here, only the charge profile is calculated via 



DFT in their work is used. More experimental data un- 
der controlled situations is needed for further testing pre- 
dictions of flow. 

We also briefly draw attention to induced-charge 
electro-osmotic flows (ICEO) @, [Hj], where the new 
model may help to explain some puzzling experimen- 
tal results (although we do not report any new simu- 
lations here). In particular, we (along with L. R. Ed- 
wards and M. S. Kilic) showed that flow reversal in AC 
electro-osmotic micro-pumps (consisting of interdigitated 
planar micro-electrode arrays) could be explained by a 
Bikerman-like model of the double layer, where the dif- 
ferential capacitance of the double layer decreases at high 
voltage . A difficulty with this interperation of the ex- 
perimental data, however, was the fact that the inferred 
ion size was far too large, whether considering a lattice 
gas or hard spheres. The problem could be alleviated by 
considering the possibility of reduction of the dielectric 
constant near the surface, and we speculated that corre- 
lation effects might further reduce the effective ion size 
in the model. From the present work, we can see that 
electrostatic correlations tend to reduce electro-osmotic 
flow while also lowering the double-layer differential ca- 
pacitance. The former effect could be wholly or partly 
misinterpreted as a sign of charge- induced thickening (i.e. 
an increase in viscosity in a highly charged double layer 
that could also reduce the net electro-osmotic flow) , while 
the latter could reduce the capacitance without invoking 
such strong steric effects with large effective ion sizes. 
Based on this evidence, it seems plausible that the new 
model might help to describe ICEO flows at high volt- 
age and high salt concentration, which have otherwise 
resisted a complete theoretical understanding |2|. 



VI. CONCLUSIONS 

We have developed a continuum model for clcctroki- 
netic phenomena that accounts for electrostatic corre- 
lations and applied this model to electro-osmotic flow 
and streaming current in aqueous electrolytes of high 
valence and high salt concentration at a flat, homoge- 
neously charged surface. The model predicts the basic 
electric double layer structure that has been observed 
in Monte Carlo simulations; namely oscillations in the 
charge density and reversal of apparent charge of a sur- 
face based on electrokinetic flow. Without any fitting 
parameters, the continuum model which also includes fi- 
nite ion size effects reproduces features of much more 
complicated theories and simulations. While the quan- 
titative agreement between this model and Monte Carlo 
or DFT simulations is only approximate, the trends are 
much closer than found with the classical mean-field thc- 
As in the case of RTIL 



ory. As in the case ot KILL [27(, it is remarkable that 
such a simple continuum theory can predict various sub- 
tle aspects of double layer structure and electrokinetic 
phenomena at the molecular scale. 
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